#######################################
######Masculinity Threat Analysis######
#######################################
###############10/10/25##################
#######################################

####Analysis Script for Gothreau & Haas, forthcoming in Journal of Experimental Political Science###
####A Replication and Extension of Willer et al. (2013) Overdoing Gender: A Test of the Masculine Overcompensation Thesis###

library(haven)
library(ggplot2)
library(dplyr)
library(broom)
library(officer)
library(flextable)
library(forcats)
library(psych)
library(survey)
library(tidyr)
library(ggplot2)
library(patchwork)

#Set Working Directory#
MTData <- read.csv("CleanMTData.csv")

#Create datasets by gender
Data_female <- subset(MTData, gender_binary == "Women")
Data_male <- subset(MTData, gender_binary == "Men")

# List of dependent variables
dvs <- c("IMMIGRATION_rescaled", "MARIJUANA_rescaled", "ELECTRIC_CAR_rescaled", "TRANS_RIGHT_rescaled", "HIRINGPOLICY_WOMEN_rescaled")

run_weighted_ttests <- function(data, label) {
  experimental_groups <- setdiff(unique(data$P_EXP), "no threat")
  results <- list()
  
  for (dv in dvs) {
    if (!dv %in% names(data) || !is.numeric(data[[dv]])) {
      warning(paste("Skipping", dv, "- missing or not numeric"))
      next
    }
    
    for (group in experimental_groups) {
      # Subset and clean data
      subset_data <- data %>%
        filter(P_EXP %in% c("no threat", group)) %>%
        select(P_EXP, WEIGHT, all_of(dv)) %>%
        na.omit()
      
      if (nrow(subset_data) < 3) next
      
      # Define survey design
      design <- svydesign(ids = ~1, weights = ~WEIGHT, data = subset_data)
      
      # Run weighted t-test
      formula <- as.formula(paste(dv, "~ P_EXP"))
      ttest <- svyttest(formula, design)
      
      # Extract test statistics
      t_val <- as.numeric(ttest$statistic)
      p_val <- ttest$p.value
      
      # Get weighted means
      group_means <- svyby(as.formula(paste("~", dv)), ~P_EXP, design, svymean, na.rm = TRUE)
      mean_no_threat <- group_means[group_means$P_EXP == "no threat", 2]
      mean_treat     <- group_means[group_means$P_EXP == group, 2]
      
      # Compute difference and estimate standard error
      diff_means <- as.numeric(mean_treat - mean_no_threat)
      se <- ifelse(!is.na(t_val) && t_val != 0, abs(diff_means / t_val), NA)
      
      # Store results
      result <- data.frame(
        gender = label,
        DV = dv,
        experimental_group = group,
        mean_no_threat = round(mean_no_threat, 2),
        mean_treat = round(mean_treat, 2),
        t = round(t_val, 2),
        p = p_val,
        std_error = round(se, 3)
      )
      
      # Significance stars
      result$sig <- ifelse(result$p < .001, "***",
                           ifelse(result$p < .01, "**",
                                  ifelse(result$p < .05, "*", "")))
      
      # Store each comparison result
      results[[paste(label, dv, group, sep = "_vs_")]] <- result
    }
  }
  
  bind_rows(results)
}

ttests_female <- run_weighted_ttests(Data_female, "female")
ttests_male   <- run_weighted_ttests(Data_male, "male")
all_ttests <- bind_rows(ttests_female, ttests_male)

####Table with New DVs (C5 in Appendix)

# Keep only comparisons where experimental group is "threat"
threat_table <- all_ttests %>%
  filter(experimental_group == "threat") %>%
  mutate(
    comparison = paste0("t = ", t, sig, ", p = ", sprintf("%.3f", p)),
    DV = factor(DV, levels = dvs)
  ) %>%
  select(gender, DV, mean_treat, mean_no_threat, comparison)

ft <- flextable(threat_table)
ft <- flextable::footnote(
  ft,
  i = 1, j = 1,
  value = as_paragraph("* p < .05, ** p < .01, *** p < .001"),
  ref_symbols = "†",
  part = "body"
)
ft <- add_footer_lines(ft, "* p < .05, ** p < .01, *** p < .001; † indicates significance in original study")
ft <- autofit(ft)

doc <- read_docx()
doc <- body_add_par(doc, "T-Test Results: Threat vs. No Threat Condition", style = "heading 1")
doc <- body_add_flextable(doc, ft)
print(doc, target = "Study3_replication.docx")

########Data Visualization####### (Figure 1 in Main Text)

dv_labels <- c(
  "IMMIGRATION_rescaled" = "Immigration",
  "MARIJUANA_rescaled" = "Marijuana",
  "ELECTRIC_CAR_rescaled" = "Electric Car",
  "TRANS_RIGHT_rescaled" = "Trans Rights",
  "HIRINGPOLICY_WOMEN_rescaled" = "Hiring Women"
)

#Means for lollipop
means_data <- all_ttests %>%
  filter(experimental_group == "threat") %>%
  select(gender, DV, mean_no_threat, mean_treat) %>%
  pivot_longer(cols = c(mean_no_threat, mean_treat),
               names_to = "condition", values_to = "mean") %>%
  mutate(condition = recode(condition,
                            "mean_no_threat" = "Control",
                            "mean_treat" = "Threat"))

#Treatment effects
effect_data <- all_ttests %>%
  filter(experimental_group == "threat") %>%
  mutate(
    treatment_effect = mean_treat - mean_no_threat,
    lower = treatment_effect - 1.96 * std_error,
    upper = treatment_effect + 1.96 * std_error,
    DV_label = dv_labels[DV]
  )

means_data$DV_label <- dv_labels[means_data$DV]

#means_data$DV <- factor(means_data$DV, levels = names(dv_labels), labels = dv_labels)
#effect_data$DV <- factor(effect_data$DV, levels = names(dv_labels), labels = dv_labels)

means_women <- means_data %>% filter(gender == "female")
effects_women <- effect_data %>% filter(gender == "female")

p1_women <- ggplot(means_women, aes(x = mean, y = fct_rev(DV_label), color = condition)) +
  geom_line(aes(group = DV), color = "gray80") +
  geom_point(size = 3) +
  scale_color_manual(values = c("Control" = "gray60", "Threat" = "darkblue")) +
  labs(title = "Women", x = "Mean Outcome", y = NULL, color = NULL) +
  theme_minimal(base_size = 20) +
  theme(legend.position = "bottom")

p2_women <- ggplot(effects_women, aes(x = treatment_effect, y = fct_rev(DV_label))) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  # 95% CI (wider)
  geom_errorbarh(aes(xmin = treatment_effect - 1.96 * std_error,
                     xmax = treatment_effect + 1.96 * std_error),
                 height = 0.25, color = "darkblue", alpha = 0.4) +
  # 90% CI (narrower)
  geom_errorbarh(aes(xmin = treatment_effect - 1.645 * std_error,
                     xmax = treatment_effect + 1.645 * std_error),
                 height = 0.1, color = "darkblue", alpha = 0.8) +
  geom_point(color = "darkblue", size = 3) +
  labs(title = NULL, x = "Treatment Effect (Threat – Control)", y = NULL) +
  theme_minimal(base_size = 20) + 
  theme(panel.grid = element_blank())

means_men <- means_data %>% filter(gender == "male")
effects_men <- effect_data %>% filter(gender == "male")

p1_men <- ggplot(means_men, aes(x = mean, y = fct_rev(DV_label), color = condition)) +
  geom_line(aes(group = DV), color = "gray80") +
  geom_point(size = 3) +
  scale_color_manual(values = c("Control" = "gray60", "Threat" = "darkblue")) +
  labs(title = "Men", x = "Mean Outcome", y = NULL, color = NULL) +
  theme_minimal(base_size = 20) +
  theme(legend.position = "bottom")

p2_men <- ggplot(effects_men, aes(x = treatment_effect, y = fct_rev(DV_label))) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  geom_errorbarh(aes(xmin = treatment_effect - 1.96 * std_error,
                     xmax = treatment_effect + 1.96 * std_error),
                 height = 0.25, color = "darkblue", alpha = 0.4) +
  geom_errorbarh(aes(xmin = treatment_effect - 1.645 * std_error,
                     xmax = treatment_effect + 1.645 * std_error),
                 height = 0.1, color = "darkblue", alpha = 0.8) +
  geom_point(color = "darkblue", size = 3) +
  labs(title = NULL, x = "Treatment Effect (Threat – Control)", y = NULL) +
  theme_minimal(base_size = 20)+
  theme(panel.grid = element_blank())

#Left = Means (Control vs. Threat), Right = Treatment Effect
(p1_women | p2_women) /
  (p1_men   | p2_men)

ggsave(path = "Visualization", filename = "Means_ATE_newdvs.jpg", width = 15, height = 11, dpi=700)
dev.off()











